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ABSTRACT 

A new integration algorithm which has the simplicity of Euler integration but exhibits 
second-order accuracy is described. In fixed-step numerical integration of differential equations 
for mechanical dynamic systems the method represents displacement and acceleration variables 
at integer step times and velocity variables at half-integer step times. Asymptotic accuracy of the 
algorithm is twice that of trapezoidal integration and ten times that of second-order Adams- 
Bashforth integration. The algorithm is also compatible with real-time inputs when used for a 
real-time simulation. It can be used to produce simulation outputs at double the integration 
frame rate, i.e., at both half-integer and integer frame times, even though it requires only one 
evaluation of state- variable derivatives per integration step. The new algorithm is shown to be 
especially effective in the simulation of lightly-damped structural modes. Both time-domain and 
frequency-domain accuracy comparisons with traditional integration methods are presented. 
Stability of the new algorithm is also examined. 

1. Introduction 

In the simulation of mechanical dynamic systems described by ordinary differential 
equations the required dynamic accuracy is often modest, especially when the real-time 
computation is utilized as part of a hardware-in-the-loop simulation. Accuracies ranging 
between 0.1 and 1 percent are considered adequate in many cases. For this reason, lower-order 
numerical integration algorithms are often employed. Also, fixed integration time steps are 
invariably used in real-time simulations in order to assure that the simulation outputs for each 
integration step occur at a fixed rate that can be synchronized with real time. In fact, the Adams- 
Bashforth second-order predictor algorithm, hereafter referred to as AB-2, is perhaps the most 
widely used method for real-time simulation. 

In this paper we consider a modified form of Euler integration which is well suited to the 
dynamic simulation of mechanical systems. It is especially effective in the simulation of 
systems with lightly-damped oscillatory modes, such as flexible structures. The method has the 
simplicity of conventional Euler integration but exhibits dynamic errors that are second order 
rather than first order in the integration step size h. Also, the dynamic error coefficients 
associated with the method are smaller than those for any other second-order method. In the 
next section we introduce the basic concept behind the modified Euler method as it is used in the 
dynamic simulation of mechanical systems. This is followed by a discussion of dynamic error 
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measures with emphasis on the frequency domain. Several example simulations are then 
introduced to demonstrate the accuracy improvement achieved when using the modifed Euler 
method instead of conventional algorithms. The stability boundaries for different versions of 
the modified Euler method are also compared with those for conventional methods. 

2. The Modified-Euler Method 

The simulation of mechanical dynamic systems normally requires the integration of a 
time-dependent acceleration A(r) to obtain a velocity V, followed by a second integration to 
obtain a displacement D. This is illustrated diagramatically in Figure 1 One method for 
implementing the required integrations is to use the forward Euler formula for the first 
integration and the backward Euler formula for the second integration. The required difference 
equations are the following: 

^n+i = + hA n i D n+ j = D n + hV n+l (1) 

Here h is the integration step size and A n , V n and D n represent the respective variables at the 
time t = nh, where n is an integer. Eq. (1) has been used in real-time simulation to achieve 
dynamic accuracy improvement over that obtained when using the forward Euler formula for 
both integrations. In Eq. (1) the first-order error associated with the forward Euler formula 
cancels the equal and opposite first-order error associated with the backward Euler formula. As 
a result the displacement D exhibits second-order accuracy with respect to the input acceleration 
A. 
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Figure 1. Paired integration to obtain velocity and displacement from acceleration. 


Both integrations become second order if we consider the velocity to be represented at a 
half-integer frame. In this case the difference equations become 

^n+l/2 = ^n-l/2 + h An , ^n+1 = &„ + hV n+l/2 (2) 

The acceleration A will, of course, usually be a function of both velocity V and displacement D, 
as well as an explicit time-dependent input U(t). In this case we can write the system state 
equations as 

V = A[D,V,U(t )] , D = V (3) 
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where in general the variables will be vectors rather than scalars. In Eq. (3) we see that the 
acceleration A n at the nth frame depends on the velocity V n at the nth frame, which is not 
available in the half-integer representation for V as utilized in Eq. (2). The best we can do is to 
employ an estimate for V n based on half-integer values V n . m-> V n-zrii etc - Then the modified 
Euler difference equations are given by 


= V^ ta + hA{D n ,V„U„) , D m = D n + hV ntm (4) 

Table 1 lists some possible candidate formulas for estimating V n . In the first formula in Table 1 
we let = Vfl-i/ 2 - This is equivalent to using conventional Euler integration rather than 
modified Euler integration for the V dependent portion of A(D ,V ,IT) , with the corresponding 
dynamic error proportional to h. The second formula for V n uses a linear extrapolation based on 
V n -\p. and Vn-3/2- It is equivalent to using AB-2 integration for the V dependent portion of A, 
with the corresponding dynamic error proportional to /i 2 . In the third formula v n is derived from 
averaging V„+i/2 and V n -\n- It is equivalent to trapezoidal integration for the V dependent 
portion of A and represents an implicit formulation, since V«+i/2 now appears on both sides of 
the left equation in (4). Later in this section we will see how this can be turned into an explicit 
formulation in many cases. Finally, the last formula in Table 1 uses a second-order predictor 
integration method to obtain $ n from V n .\/2 and the derivatives V n .\ and V n - 2- It produces a 
local truncation error in proportional to h? and therefore permits the full accuracy of the 
modified Euler method to be realized. 


A A 

Table 1. Methods for Estimating V n in A(D n ,V n ,U n ) 


Euler 

A 

Vn = V n . w 

AB-2 

V = —V - —V 
n 2 nAi 2 2 ' l 3/2 

Trapezoidal 

A ^n+l/2 + V n -ifl 

n " 2 

Predictor 

Integrator 

V = V + hi— V -—V ) 

v n v n-\a "tg •'n., g %- 2 ' 


Before considering some examples of the application of the modified-Euler methods 
described here, we consider some dynamic error measures for examining comparitive accuracy 
of different integration algorithms. 
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3. Integrator Error Measures 


Consider the solution of the state equation dyldt =/(r) using a numerical integration formula 
forjn+i in terms of y n and the derivative/. Furthermore, lety[(n+l)/i] and y[nh] represent the 
exact solution of the continuous system at the times t = ( n+\)h and nh, respectively. Then we 
can then write 

y n+ i - y n - y - y m ■ e ifi k)fik+l 


Here the term -e f fyph k+ 1 represents the local truncation error associated with the integration 
method of order k and /(*) is the kth time derivative of/at t = nh [1]. For example, k = 1 and e { 
= 1/2 for Euler integration; for AB-2 integration k = 2 and = 5/12. We now take the Z 
transform of Eq. (5) and divide by z - 1 to obtain 

e, h k+l F ( *>* (z) 

7*(z) = Y f *{z) - 7 ( 6 ) 


1 


Here Y re j*{z) is the Z transform of the exact solution, y[nh]. Next we consider the case of 
sinusoidal data sequences by replacing z with ei ah . We also note that F(*)*(e/® A ) = 
( j(o) k F*(ei° )h ), i.e., the Fourier transform of the kth derivative of a function is equal to the 
Fourier transform of the function multiplied by (jco) k . After dividing the resulting expression by 
F*, we have 

TV'**) = y re *(« /fl *> _ «/*(/«*>* (?) 

F*(e iah ) F*{e ^ mh ) e ’ ah - 1 

The term Y*/F* is simply the sinusoidal transfer function, ///(e/^), of the numerical integrator. 
The term Y re f/F* = \/jco, the transfer function of an ideal integrator. If we now approximate 
d ah -\ by jcoh, Eq. (7) becomes the following: 



1 - ejijahf 

jco 


t- , coh« 1 

j(0[ 1 + ej(jcoh) ] 


( 8 ) 


Here ej is the integrator error coefficient and k is the algorithm order. To illustrate the application 
of our integrator transfer function model, we consider the simulation of a linearized dynamic 
system with transfer function H(s). For the case of sinusoidal inputs of frequency ft), the 
transfer function becomes H(jco). When the continuous system is simulated with a single-pass 
integration method, the sinusoidal transfer function of the digital simulation is simply given by 
H(\/Hj*), where H { * is the transfer function of the digital integrator. For coh « 1, l/Hj* can be 
approximated by jco[ 1 + e } (jcoh)^ in accordance with the integrator model of Eq. (8). Thus the 
formula for the transfer function of the digital system in simulating the linear system with transfer 
function H(s) is given by 


H V**) = = H{M\+ej(j(oh ) k ]} , coh« 1 (9) 
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For example, consider a first-order linear system with eigenvalue A and transfer function given 
by H(s) = 1/(5 - A). The transfer function for sinusoidal inputs becomes 

H(jco) = -r-L- (10) 

jco - A 


Then the digital system transfer function for sinusoidal input data sequences is given approxi- 
mately by 

H*^ 0 *) = — r , jco«l (11) 

ja>[l + e { (jcoh) ] - A 


We note that the characteristic root (eigenvalue) A for the continuous system is given by the value 
of jo in Eq. (10) which makes the denominator vanish. It follows that the equivalent 
characteristic root A* for the digital system is given by the value of jco which makes the 
denominator of Eq. 11) vanish. Replacing jco by A* in Eq. (11) and setting the denominator 
equal to zero, we can write 

A* = A - X*e I (X*h) k 

For | A/il « 1, A* = A to order /i k . Then we can replace X* by A on the right side of the equation 
and obtain 

e. = = . e,(Xh) k , IA/tl « 1 (12) 

K A 1 


Here e 1 represents the fractional error in the digital system characteristic root. We recall that the 
transfer function for any finite order linear system with distinct roots can be represented as the 
sum of first-order transfer functions of the form l/(s - A), where the charcteristic roots may be 
real or complex. It follows that Eq. (12) can be used to estimate the error in each characteristic 
root in the digital system simulation of any order linear system. 

From the digital transfer function formula in Eq. 11) we can write 


from which 


tf*(e /a *) 


1 


(jco- A) 



ejjco(jcoh) k 
ja>-X . 


Hjjco) 

e.jco(jcoh) k 

1 + — : ; — 

jco- X 


wVg) - Hm _ yoKW* _ 0>h<<l 

H(jco) jco- X 


(13) 


Here ( H * - H)/H represents the fractional error in digital system transfer function. For coh « 1 
it is evident that this fractional error will be small in magnitude compared with unity. In this case 
it can be shown that the real part of (H* - H)!H is approximately equal to the fractional gain error 
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of the digital transfer function and the imaginary part is approximately equal to the phase error 
[2]. We note that the transfer function for any finite-order linear system can be written as the 
product of individual pole and zero factors of the form ( s - X), where again X can be either real or 
complex. It is then straightforward to show that the fractional error in the overall digital transfer 
function is approximately the sum of the individual errors given by Eq. (13) for each factor [2]. 
It follows that both gain and phase errors of the overall digital system transfer function for 
sinusoidal inputs are proportional to efj(ah) k . 

Thus for single-pass integration methods Eq. (12) and (13) represent simple approximate 
formulas for both characteristic root and transfer function errors. For a given integration 
algorithm the errors are directly proportional to the integrator error coefficient e t for that 
algorithm. Table 2 lists and k for the algorithms considered in this paper, including the 
modified-Euler method, which has the smallest error coefficient {e f = 1/24). 


Table 2. Error Coefficients for Integration Methods 

Integrator transfer function = J7 / *(e , °*) = Y , coh« 1 

jo)[ 1 + e^jcoh) ] 


Euler 

1 

7 

1 

AB-2 

5 

17 

2 

Trapezoidal 

1 

'17 

2 

Modified Euler 

1 

24 

2 


All of the above algorithms in Table 2 are explicit except trapezoidal, which is implicit. 
An explicit method with the same asymptotic accuracy can, however, be realized with the two- 
pass Adams-Moulton (AM-2) algorithm. In this method the first pass employs AB-2 integration 
to obtain an estimate $ n+1 of the next state. This is then used in the trapezoidal formula to 
compute the corrected y n +\ . The local truncation error associated with y n+1 is of order h}, 
which ensures that the asymptotic accuracy of order H 2 for the corrected y n+ \ will be the same as 
that for implicit trapezoidal integration. 

4. Specific Examples 

We now turn to some specific examples to compare the accuracy of modified Euler 
integration with traditional algorithms of second order. We consider first a simple linear 
dependence of the acceleration A in Eq. (3) on the displacement D and velocity V. This leads to 
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the following state equations, which for convenience have been written in terms of the undamped 
natural frequency (On and damping ration £ of the second-order system. 

V = ccfc(U n -D n ) - 2£cq,V , t> = V (14) 

From Eqs. (2) and (4), with the trapezoidal formula from Table 1 used for fyt+i, we obtain the 
following difference equations for the modified Euler formulation: 

v n+m = V n _ m + o? n h{U n -D n ) - Cco n (y n+m + V nAa ) , D n+l = D n + hV^ m 


After solving the first equation for V„+i/2, we have the following explicit equations. 


where 


V„ m = c,v„ B + CJU.-D.) , 

D mi -D„ + hv „, 0 

(15) 

r - c = 

(£h 

(16) 

1 1 + C<H,A ’ 2 

i + £<**,* 



Here the constants C\ and C2 can be precomputed. From Eq. (15) it follows that the ongoing 
simulation run then requires only 3 adds and 3 multiples per integration step. 

Figure 2 shows plots of the solution error when using modified Euler integration to 
compute the response of a second-order system with £ = 0.25 to a unit step input. The initial 
conditions are given by *(0) = y(0) = 0. Shown in the figure are error plots for three of the 
damping methods listed in Table 1, including the trapezoidal damping used to derive Eqs. (15) 
and (16). For comparison Figure 2 also shows the step-response errors when AB-2 integration 
is used. In all cases the step size is given by (o n h = 0.25. The startup problem associated with 
AB-2 integration (the initial states at t = -h are not specified) is solved by using Euler integration 
for the first step. In the case of modified Euler integration the first step which computes y m 
from yo uses a step equal to hi 2. The figure clearly shows the superior accuracy ofthe modified 
Euler method, with the scheme using second-order predictor integration to estimate y n producing 
the smallest errors. 

The second example considered in this section is the simulation of the full nonlinear flight 
equations of an aircraft. Since the largest characteristic roots for the rigid airframe are nm y 
those associated with the short-period pitching motion, we will only consider symmetric flight, 
i e the longitudinal equations of motion, in our example simulation. The conclusions regarding 
dynamic errors can be safely extrapolated to the full six-degree-of-freedom case. For this 
simulation the translational equations of motion are written with respect to flight-path axes, while 
the rotational equations of motion are written with respect to body axes [3]. Then the velocity 
state variables become total aircraft velocity V p , angle of attack a, and pitch rate Q. e 
displacement state variables are altitude H, pitch angle <9, and horizontal distance X. The 

velocity state equations are given by 

a = Q+Ltil- . O = — (17) 
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v p m 
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Figure 2. Unit step response errors in simulating second-order system, £ = 0.25, (O n h = 0.25. 
and the displacement state equations by 

0 = Q , H = V p sin(@-a) , X = V p co%{0- a) (18) 

Here F wx and F wz are the external force components along the x and z flight-path axes, 
respectively, and M is the moment about the y body axis; m and lyy represent, respectively, the 
aircraft mass and pitch-axis moment of inertia. The following formulas were used to represent 


the external forces and moment: 

F wx =-qS (,C Do +C Dc i Cl ) - gsin(0-a)+^-cosa (19) 

F wz = qS (C L +C Lg 5 e ) + gcos(0-a)-£sina (20) 

M = qcS (C Mq + C Ma a + c M q -^Q + C Md 2 vT ® + C M s fi ) 

where 

q = dynamic pressure = (22) 

and 

C L = lift coefficient = C^+C^a (23) 
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In these equations S is the aircraft wing area, g is the gravity acceleration, T is powerplant thrust, 
5 e is elevator displacement, and c is the mean aerodynamic chord. The various C s represent 
aerodynamic coefficients and stability derivatives in accordance with the subscripts. In a full 
flight-envelope simulation these will be nonlinear functions of other variables such as V p 
(through Mach number dependence), a, S e , and h. 

Based on the way in which modified-Euler integration was introduced in Section 2, the 
velocity states V p , a and Q in Eq. (17) would be represented at half-integer frames, with the 
position states 0 , H and X represented at integer frames. For the nth integration frame this 
results in the computation of the n+ 1/2 velocity state from the n- 1/2 velocity state, followed by 
computation of the n+1 position state from the n position state using the n+ 1/2 velocity state just 
obtained. However, from Eq. (17) it is apparent that it would be better to represent the angle of 
attack a at integrer frames, even though it is derived from a velocity state equation. This is 
because the dominant term on the right side of Eq. (17) affecting the high-speed dynamics is the 
pitch-rate Q, which is represented at half-integer frames. The other term in Eq. (24), F wz JmV p , 
is the negative of the flight-path-axis pitch rate, and is generally much smaller in magnitude than 
Q. For this reason we have chosen to represent a at integer frames in the modified Euler 
mechanization of the flight equations. Since the force (and hence acceleration) term F wz is 
computed and therefore represented at integer frames, it is necessary to compute an estimate of 
F wz at the n+l/2 frame in the modified Euler integration of (X n to obtain a n + 1 . This is easily 
accomplished using the first-order extrapolation formula F WZ/J+1/2 = (3/2 )F WZn - (1/2 )F WZn . v The 
actual difference equations used to solve (17) through (23) with modified Euler integration are 
presented in a previous paper by the author [4]. 

As a specific example we consider a business jet flying at 40,000 feet at a speed of Mach 
0.7 [5]. For the above flight condition the undamped natural frequency of the short-period mode 
is about 3 rad/sec and the damping ratio is 0.4. In the example simulation we let the step size h = 
0.1 second. This makes a>nh = 0.3 for the short-period motion, which should yield the moderate 
accuracies normally associated with a real-time simulation. We consider the aircraft response to 
the input function shown in Figure 3, which is a step elevator displacement with a one second 
rise time. Use of this input function tends to reduce the large transient errors caused by step 
inputs when predictor integration algorithms are used. It is also probably more typical of an 
actual transient input. The simulation is started at t = 0 with the aircraft in level equilibrium 
flight. In order to make the example more representative of an ongoing simulation, the step input 



Figure 3. Delayed, finite rise-time step input. 
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is delayed for 0.3 seconds (three integration steps for h = 0.1) after the initial time t - 0. Figure 
4 shows the error in pitch angle versus time for AB-2 integration and for modified Euler 
integration using the predictor integration of Table 1 for the integer velocity estimate. We note 
that the modified Euler method is an order of magnitude more accurate. 
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Figure 4. Aircraft pitch angle error for the input function of Figure 3; h = 0.1 seconds. 


5. Stability Considerations 

In addition to considering the dynamic accuracy associated with different numerical 
integration methods, it is important to consider the stability of the methods. This is usually done 
by considering the stability boundary in the complex XJt plane. These boundaries are shown in 
Figure 5 for modified Euler integration used to solve Eq. (14) with the various methods for 
computing the velocity estimate, V^, as presented in Table 1. Also shown in Figure 5 is the 
stability boundary for the AB-2 predictor method, as well as that for the AM-2 predictor- 
corrector method. In the latter case the stability region has been reduced by a factor of two to 
take into account that AM-2 is a two-pass method. Any values of Xh lying outside the boundary 
shown for a given method (the boundaries are symmetric with respect to the real axis) will lead to 
instability. From the figure it is evident that the modified Euler method with trapezoidal 
integration for the damping term exhibits the largest stability boundary. Note also that the 
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stability boundary for all of the modified Euler methods lies on the imaginary axis. This means 
that modified Euler integration, when used to simulate systems with pure imaginary roots, as in 
the case of undamped oscillatory modes, will also exhibit pure imaginary roots corresponding to 
zero damping. This is true regardless of the integration step size h and is the reason why the 
modified Euler method is especially effective in simulating lightly-damped dynamic systems. 



Figure 5. Stability boundaries for modified Euler and other second-order integration methods. 


In the second-order system example considered in Section 4 we were able to use 
trapezoidal integration for the damping term in the modified Euler mechanization because the 
damping was linear. This in turn permitted us to construct an explicit, single-step formulation 
represented by Eqs. (15) and (16). When the the dependence of acceleration on velocity is 
nonlinear, this is no longer possible. Yet it would be advantageous for stability reasons to still 
use a trapezoidal implementation. 

The nonlinear dependence of the acceleration A on the velocity V in Eq. (3) can often be 
expressed in terms of VBA/BV, where BA/BV is not a function of V, or at worst is only slightly 
dependent on V. For example if A represents dQIdt, the time derivative of pitch rate Q in the 
flight equations, then BA/BQ is proportional to the aerodynamic stability dericvative Cmq, i-e., 
the dimensionless pitching moment due to dimensionless pitch rate. Cmq is normally 
independent of Q, although it may be dependent on other variables such as Mach number. Also, 
the overall BA/BQ in this case will be independent of Q. Letting V be a scalar which represents 
the angular velocity Q, we can rewrite Eq. (14) as follows: 


710 


ORIGINAL PAGE IS 
OF POOR QUALITY 








(24) 


V = C Q [D,U(t)]+ C x [D,U(tW 

where Co + C\V = A and Ci = dA/dV. Now, when mechanizing the modified-Euler difference 
equations (4) we can compute ty n , the estimate of V at the nth frame, by the formula 

y - l(y +y ) (25) 

From Eqs. (24) and (25) the difference equation for V n +i /2 in Eq. (4) then becomes 

F„.,„ = V„ m + «C 0 (D„,t/„) + (26) 

With respect to the velocity state V this equation clearly represents implicit trapezoidal 
integration. However it can be solved to obtain the following explicit formula for V n+ i/2'- 


^n+l/2 ~ 


(1+^/2)^ + AC 0 
I-ZiCj/2 


This formulation, i.e., the use of trapezoidal integration for the damping term, expands very 
substantially the stability region in the Xh plane compared with the use of the predictor formula 
for ^ n , as we have seen in Figure 5. It can also reduce appreciably the dynamic errors 
following transient inputs. The extra required computation is modest and consists mainly of an 
additional division. 

In deriving Eq. (27) we have assumed that V is a scalar, whereas V will in general be a 
vector. In this case dA/dV will be a matrix, which must be inverted to obtain the explicit formula 
for Vn+i/ 2 - Fortunately, the critical terms in this matrix in the case of the flight equations are the 
diagonal terms, in which case simple formulas similar to Eq. (27) involving only the diagonal 
terms can be derived. In the longitudinal flight equations (17) through (23), for example, an 
equation similar to (27) can be written for Q n + 1 / 2 , where C\ is proportional to the stability 
derivative Cmq- 


6. Conclusions 

We have shown that mechanical dynamic systems are well suited to a modified Euler 
integration method which computes displacement and acceleration variables at integer frame times 
and velocity variables at half-integer frame times. Examination of asymptotic formulas for 
characteristic root and transfer function errors associated with a linearized version of any 
nonlinear mechanics problem shows that the modified Euler method is at least twice as accurate 
as any other known second-order algorithm. For the usual case where the acceleration is a 
function of velocity, there are a number of candidate methods for computing the required velocity 
estimates at integer frames from the velocity as computed at half-integer frames. A second-order 
predictor integration formula produces the most accurate integer-frame velocity estimate; an 
estimate based on the equivalent of trapezoidal integration produces the most stable simulation. 
Neither estimate requires any additional derivative evaluations, and the predictor formula can be 
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used to produce output displacements at half-integer as well as integer frame times in a real-time 
simulation, i.e., at double the integration frame rate. The modified Euler method is particularly 
effective in simulating systems with lightly damped modes, since modes with zero damping in a 
continuous system generate modes with zero damping in the modified Euler mechanization, 
regardless of the integration step size. The modified Euler method also has a simple and accurate 
startup procedure and is completely compatible with real-time inputs. Two examples, a second- 
order linear system and a sixth-order nonlinear flight simulation, have been used to demonstrate 
the superior accuracy of the modified Euler method. 
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